1.从FNA文件中提取指定区间的序列

方法1：使用seqkit工具（推荐）

# 安装seqkit (如果尚未安装)
conda install -c bioconda seqkit

# 提取指定区间 (例如：从100到200的位置)
seqkit subseq --chr "序列ID" -r 100:200 input.fna > output.fna

方法2：使用samtools faidx

# 首先建立索引
samtools faidx input.fna

# 然后提取区间
samtools faidx input.fna "序列ID:100-200" > output.fna

方法3：使用Python脚本

from Bio import SeqIO

def extract_sequence(fna_file, seq_id, start, end, output_file):
    for record in SeqIO.parse(fna_file, "fasta"):
        if record.id == seq_id:
            subseq = record.seq[start-1:end]  # Python是0-based，生物坐标是1-based
            with open(output_file, "w") as f:
                f.write(f">{record.id}_{start}_{end}\n{subseq}\n")
            return
    print(f"序列ID {seq_id} 未找到")

# 使用示例
extract_sequence("input.fna", "chr1", 100, 200, "output.fna")

方法4：使用bedtools

# 首先创建一个BED文件定义区间
echo -e "序列ID\t99\t200" > region.bed  # BED是0-based半开区间

# 然后提取序列
bedtools getfasta -fi input.fna -bed region.bed -fo output.fna

注意事项

1. 区间定义：不同工具对区间定义可能不同(0-based或1-based，闭区间或半开区间)
2. 序列ID：确保输入的序列ID与文件中的完全匹配
3. 大文件处理：对于大基因组文件，建议先建立索引(如samtools faidx)
4. 反向互补链：如果需要提取反向互补链，可以添加参数如--reverse-complement(seqkit)或处理Python中的序列

2.从FASTA文件 map_gene.fa 中保留最长的200条序列

方法2：使用 seqkit（更高效，推荐）

seqkit sort --by-length --reverse map_gene.fa \
  | seqkit head -n 200 > top200.fa

步骤解释：

1. seqkit sort 按序列长度降序排序。
2. seqkit head 提取前200条记录。

安装 seqkit：

wget https://github.com/shenwei356/seqkit/releases/download/v2.3.0/seqkit_linux_amd64.tar.gz
tar -zxvf seqkit_linux_amd64.tar.gz
sudo mv seqkit /usr/local/bin/

方法3：纯 awk + 基础命令

awk '/^>/ { 
        if (header) print header ORS seq
        header=$0
        seq=""
        next
     } 
     { seq = seq $0 } 
     END { if (header) print header ORS seq }' map_gene.fa \
  | awk 'BEGIN {RS=">"; FS="\n"} NR>1 { 
        len=0; seq=""
        for (i=2; i<=NF; i++) { len+=length($i); seq=seq $i }
        print len, ">"$1, seq
     }' \
  | sort -k1,1nr \
  | head -n 200 \
  | cut -d' ' -f2- \
  | tr ' ' '\n' \
  | awk '/^>/ {print; next} {gsub(/.{80}/,"&\n"); printf "%s", $0}' > top200.fa

步骤解释：

1. 合并多行序列为单行。
2. 计算每条序列的长度并附加到行首。
3. 按长度降序排序，取前200条。
4. 恢复FASTA格式，每80字符换行。

3.批量修改FASTA文件序列ID的几种方法

方法1：使用Linux命令 awk

awk 'BEGIN{FS=OFS="\t"; while(getline < "chrom_map.txt") {map[$1]=$2}} 
     /^>/ {id=substr($0,2); split(id,a," "); old=a[1]; 
          if (old in map) { $0 = ">"map[old] substr(id,length(old)+1) } } 
     1' input.fna > output.fna

说明：

1. BEGIN块读取chrom_map.txt到字典map。
2. 处理以>开头的行，提取旧ID并替换为新ID，保留后续描述。
3. 1表示打印所有行。

方法2：使用生物软件SeqKit（推荐）

seqkit replace -p "^(\S+)" -r '{kv}' -k chrom_map.txt input.fna -o output.fna

步骤：

1. 安装SeqKit（若未安装）：

   conda install -c bioconda seqkit  # 通过conda
   # 或
   brew install seqkit               # 通过Homebrew (macOS)
   

2. 运行替换命令，-p匹配ID，-k指定映射文件。

方法3：使用Python脚本

import sys

def main():
    chrom_map = {}
    with open('chrom_map.txt', 'r') as f:
        for line in f:
            old, new = line.strip().split('\t')
            chrom_map[old] = new

    with open('input.fna', 'r') as infile, open('output.fna', 'w') as outfile:
        for line in infile:
            if line.startswith('>'):
                parts = line[1:].split(maxsplit=1)  # Split into ID and optional description
                old_id = parts[0]
                new_id = chrom_map.get(old_id, old_id)
                new_line = f'>{new_id}'
                if len(parts) > 1:
                    new_line += ' ' + parts[1]
                outfile.write(new_line + '\n')
            else:
                outfile.write(line)

if __name__ == '__main__':
    main()

运行：

python replace_ids.py

说明：

1. 读取映射文件到字典。
2. 处理FASTA文件，保留ID后的描述信息。